Impact of host phonons on interstitial diffusion

The net effect of host phonons on interstitial diffusion has remained as a fundamental knowledge gap in our current theories since the motions of the host atoms and interstitials were coupled in these theories. Here we study this effect through molecular dynamics simulations of hydrogen diffusion in palladium, in which the motions can be decoupled through pinning the host atoms. Mathematically this decoupling corresponds to expanding the total diffusion coefficient into a Taylor series, which separates the phonon contribution from the intrinsic interstitial jumping. Our results clearly show that palladium phonons significantly promote hydrogen diffusion. The phonon contribution, being linear with temperature at high temperatures and exponential at low temperatures, is fitted with Brownian motion model. The total diffusion of interstitials can be understood as the intrinsic interstitial jumping in a pinned host plus phonon-induced Brownian diffusion. The generality of our findings is validated by examining the motion of lithium in manganese oxide and carbon in iron.


Impact of host phonons on interstitial diffusion
Chunguang Tang 1,2* , Gang Sun 3 & Yun Liu 1,2 The net effect of host phonons on interstitial diffusion has remained as a fundamental knowledge gap in our current theories since the motions of the host atoms and interstitials were coupled in these theories. Here we study this effect through molecular dynamics simulations of hydrogen diffusion in palladium, in which the motions can be decoupled through pinning the host atoms. Mathematically this decoupling corresponds to expanding the total diffusion coefficient into a Taylor series, which separates the phonon contribution from the intrinsic interstitial jumping. Our results clearly show that palladium phonons significantly promote hydrogen diffusion. The phonon contribution, being linear with temperature at high temperatures and exponential at low temperatures, is fitted with Brownian motion model. The total diffusion of interstitials can be understood as the intrinsic interstitial jumping in a pinned host plus phonon-induced Brownian diffusion. The generality of our findings is validated by examining the motion of lithium in manganese oxide and carbon in iron.
Lattice vibrations or phonons play a major role in many physical properties, such as thermal and electrical conductivity, of condensed matter systems. Their effect on mass transport, i.e., atomic diffusion, however, is not well understood despite studies for decades. It is well known that the diffusion in solids occurs via atomic jumping and follows the Arrhenius equation. In particular, based on the established random walk theory 1,2 the diffusion coefficient for interstitials can be written as where E a is the activation enthalpy for the jump, k B is the Boltzmann constant, and D 0 is a constant virtually independent of temperature T and combines the effects of solid structures, interstitial vibration frequency, and an entropy change term associated with the jump. In the early versions 2,3 of the random walk theory, the effect of host phonons was not considered 4 . Later, a number of classical models 4-7 replaced the vibration frequency of the diffusing atoms with an effective frequency which involves the many body effect from the lattice phonons. The phonon effects on atomic diffusion were also discussed 8,9 based on the quantum theory 10 . In all these theoretical studies, however, the motions of the host atoms and the jumping atom are coupled and hence the net effect or importance of phonons is not clear. In other words, although one can compute the diffusion coefficients with lattice phonons included based on these theories 11 , it is not easy to know qualitatively whether lattice vibrations positively or negatively contribute to the diffusion and quantitatively how significant the contribution is. This represents a knowledge gap for our fundamental understanding of solid diffusion, although the random walk theory has been established for decades.
The advance of our computational capacity and relevant algorithms allows us to study the above problem via simulations. Here we address the effect of lattice phonons on interstitial diffusion, using hydrogen (H) in facecentered-cubic (fcc) and amorphous palladium (Pd) as an example, based on molecular dynamics simulations. We consider the amorphous phase here because of its very different structural orders as compared to traditional crystalline solids. Hydrogen in solids has a long history in the interest of scientists 12,13 , and metal hydrides for hydrogen storage have seen extensive research activities over the past decades 14 . Pd is one of the first metals of which the capacity for H absorption were discovered 15 and is still widely applied in the hydrogen energy sector 16,17 , for example, as effective hydrogenation catalysts, hydrogen purification filters and hydrogen storage media. We validate the generality of our findings by examining two technically important interstitial systems, namely, carbon (C) in iron (Fe) and lithium (Li) in manganese oxide ( MnO 2 ). www.nature.com/scientificreports/ To investigate the amorphous phase, a Pd-H system with 4000 Pd atoms and 200 H atoms was well liquefied at 2100 K and quenched to 300 K at 10 14 K/s, and the resulting structure was confirmed to be amorphous based on its Pd-Pd pair distribution function. For the crystalline system, we randomly put 200 H atoms at the octahedral interstitial sites of an fcc structure of 4000 Pd atoms. We relaxed the obtained crystalline structures at various temperatures to obtain the equilibrium system volume. These computations were carried out using the NPT (constant atom number, pressure, and temperature) ensemble by setting the normal pressures to zero. Using the NVT (constant atom number, volume, and temperature) ensemble at various temperatures with corresponding system volume obtained from the NPT ensemble, we then computed the mean squared displacements (MSDs) of H and Pd atoms in the obtained structures for various times after an initial relaxation for 50 ps or so. The MSD sampling time ranges from 0.1 to 200 ns, depending on temperature, Pd phase, and whether Pd atoms were pinned. The varying sampling time is to ensure the motion of H in a given situation becomes stable. Each MSD data point in this work represents the average over about ten independent computations. More details of MSD sampling can be found elsewhere 19 . The MSD data were used to compute the diffusion coefficient based on the Einstein relation where t is the time and r 2 (t) is the ensemble average of MSDs of diffusing atoms. For H diffusion, only MSD data higher than 10 Å 2 were used to make sure the system is in the diffusive regime. Pd diffusion is much slower, and hence we only sampled the MSD at a fixed time length (100 ps) for implicating its mobility. The computations for carbon in crystalline iron are similar to the above. The atomic interactions of the Pd-H system were described by an embedded-atom-method-based potential 20 . Simulations 21 based on this potential has confirmed the concentration dependence of hydrogen diffusion in metals as observed experimentally. All the simulations in this work were performed using the code LAMMPS (ref. 22 ) with Nosé-Hoover thermostat and barostat, and the timestep was set as 1 fs.
We started our study with the thermodynamics of Pd matrix. For fcc Pd, its simulated melting point was found to be ∼1100 K. For the amorphous phase below the glass transition temperature, T g , Pd atoms essentially only vibrate about their local favorable positions within the time scale of interest to this study (Fig. 1a). Above T g , the MSD of Pd increases monotonically with temperature. During the simulations, crystallization of Pd at elevated temperatures was avoided by selecting a reasonably short time scale (but long enough for observing stable H diffusion). We note the existence of a critical temperature ( T c ) in the supercooled region, as also reported for Pd in multiple-component amorphous alloys 23 , which corresponds to the turning point of the Arrhenius curve and signifies the change in the diffusion mode from liquid-like motion to solid-like hopping upon cooling 23 . For the purpose of this work, we limit our discussion on H diffusion in the amorphous phase below T c .
We study the net effect of host phonons by comparing H diffusion with the vibrations of Pd atoms enabled and disabled, respectively. While impossible in experiments, this strategy is convenient in simulations and useful for tackling complex problems 24 . We first simulated the diffusion of H atoms with the Pd atoms fixed (by setting the force on them to zero) after the system was annealed for some time. Note that in this case, the kinetic energy of H is still connected to an external thermostat in simulations. From a theoretical perspective, in this case the collision between H and Pd is elastic since the collision does not transfer kinetic energy to Pd. As indicated by the open circles in Fig. 1b-c, H diffusion in this case fits well into the Arrhenius equation. As can be seen, like the random walk theory, our simulations based on the static Pd matrices predict a well-behaved Arrhenius diffusion of H. We note that in the random walk theory lattice relaxation during the jumping of interstitials is possible. Here the host matrix is relaxed for some random time before it is fixed.
Next we computed D of H without pinning the Pd atoms. As indicated by the solid circles in Fig. 1b-c, H diffusion in this scenario also follows the Arrhenius relationship. Figure 1c also shows the experimental 18 D of hydrogen from ∼ 530 to ∼ 910 K, which differs from our calculations within one order of magnitude. Other measurements 25,26 reported the values of D near 300 K to be ∼ 1.6 − 1.9 × 10 −7 cm 2 /s, which are close to the extrapolation of the previous measurement 18 and differ from our calculations by more than one order of magnitude. Such a difference mainly results from the higher activation energy reported by experiments, which is not surprising in view of the trapping of hydrogen by lattice defects in experiments 27 . The accuracy of the EAM potential used in this work may also contribute to the difference, but overall the calculations resemble the trend of experiments. Figure 1b-c clearly shows a general and significant difference in the diffusion coefficients for fixed and mobile Pd. We label the diffusion coefficients for the case of fixed Pd as D j since here only H jumping effect contributes to the diffusion. Figure 1d shows that the ratio D j /D decreases as temperature decreases and reaches below 30% at room temperature. To investigate the effect of Pd phonons on H diffusion, we split the total diffusion coefficient D as where D j by itself has the Arrhenius form D j = D j 0 exp(−E j a /RT) , and D ph represents the contribution of Pd phonons. The philosophy behind this treatment is as follows. The Pd phonon effect is related to the mass of Pd ( m Pd ). Imagine H diffuses in a hypothetical Pd isotope with infinite mass ( w = 1/m Pd = 0 ). In this case, the vibration amplitude of Pd becomes zero and D reduces to D j . As w increases, the phonon effect emerges and contributes to H diffusion. This led us to formally write D as a Taylor series In the following we study D ph from the perspective of Brownian motion. Brownian motion is generally regarded as a fluid phenomenon, but similar phenomena, such as thermal noise in electric conductors 28 , exist in solids and all of them follow the fluctuation-dissipation theorem 29 . The splitting of D may raise concerns since the Brownian motion model is space/time continuous while the interstitials move in a space-discrete force field determined by the host structure. We note equation 3 implicitly treats the force field created by the host atoms as where ε f corresponds to the host atoms at their pinned positions and the fluctuating augmentation ε v arises from host phonons. The ε f field results in a discrete diffusion component ( D j ) and the Brownian contribution comes only from ε v which is continuous in space/time. Such a treatment is indeed consistent with the solid-liquid phase transition: as temperature increases, the ε v component increases and eventually the force field becomes liquid-like.
The diffusion coefficient of Brownian particles suspended in a fluid can be written as (see reference 30,31 or the Supplementary Information) where E k and a are the average kinetic energy and size of the Brownian particles, respectively, and µ is the liquid viscosity. A more familiar form for equation (6) is www.nature.com/scientificreports/ equipartition theorem. In view that H atoms are in thermal equilibrium with Pd phonons, we can directly apply Equation (6) to D ph if we know the average energy of Pd phonons E ph . In contrast to liquid of which µ is sensitive to temperature, here we can assume constant µ for Pd phonons. In other words, we assume a constant friction coefficient for the motion of H in the phonon fluid. At low temperatures 32,33 , E ph is approximately proportional to T n , where n is some constant. Substituting this into E k = 3E ph /2 , we can write At high temperatures, the equipartition theorem gives E ph ≈ k B T . Note the actual E ph (= T 0 CdT ) is smaller than k B T by some roughly constant number because the heat capacity per phonon, C, at low temperatures is lower than its high temperature limit k B . In this case, we can write the phonon contribution as where c is a negative constant dependent on the material.
As shown in Fig. 2a, at high temperatures D ph of both the amorphous and the fcc phases follows a linear behaviour with a negative intercept at T = 0 . We note that for the amorphous phase the linearity holds for a range of temperatures above T c , which we attribute to the cancellation effect as at T > T c both D j and D deviate from the Arrhenius relationship for solids. At low temperatures, a function proportional to T n (n = 4 for amorphous and 5 for fcc) seems to fit D ph well. The fitting parameters n = 4 and 5 are based on the Debye approximation 32 and the finite temperature field theory 33 , respectively. However, as shown in the log-scale inset of Fig. 2a, a constant n actually does not fit well with the data, which contrasts with the established models. To build a theory for this observation is beyond the scope of this work, but we note the results are logically not surprising since, if n can increase from 1 at high temperatures to 4 or 5 at low temperatures, it is reasonable to assume higher n at even lower temperatures. Inspired by this, we found the D ph data at low temperatures fit well with an exponential function, as shown in Fig. 2b. We note that D ph = D − D j can be approximated as D as T → 0 because, with a higher activation energy, D j is an infinitesimal of higher order than D. This explains the exponential fitting of D ph from the mathematical perspective.
For some comments to the above results, firstly we ignored the zero point vibrations in the above derivations, which only adds a small constant D ph that is important as T approaches zero. Secondly, the quantum tunneling effect for H diffusion 10 , which is important for temperatures below 200 K 34 , was not included in this work. Thirdly, we expect the phonon contribution is a general effect and hence also affect other interstitials. Indeed, a recent molecular dynamics study showed that the permeability of He atoms and H 2 molecules in amorphous silica decreases if the thermal motion of silica is forbidden 35 . To extend our findings to larger atoms, we studied the movement of carbon in fcc iron at 1500 K and that of Li + ions in Li 0.5 Mn 2 O 4 at 1000 K, and the results (in Supplementary Information) confirmed the impacts of host phonons on interstitial motion. Finally, although atom-pinning decouples the motions of the host and interstitial atoms, it is a theoretical approach without an experimental counterpart. Hence, we also examined the phonon effects by only changing the mass of Pd atoms (without pinning Pd). This is equivalent to using Pd isotopes in experiments although in simulations one can arbitrarily set the mass of Pd. We used the fcc phase at 500 K as an example. Figure 3 shows that the diffusion coefficient of H decreases as the mass of Pd increases (or, equivalently, the phonon energy decreases), which supports our findings based on Pd-pinning.
In summary, based on molecular dynamics simulations we studied the effects of host phonons on interstitial diffusion, using hydrogen (H) within face centered cubic and amorphous palladium (Pd) as an example. Compared previous theoretical studies that coupled the phonons with the interstitial motions, this work decouples the two by pinning the host atoms and hence clearly reveals the net effects of host phonons. It was found that Pd phonons significantly promote H diffusion by causing Brownian-like diffusion of H atoms and dominate H diffusion below T m /2 ( T m being the melting point of Pd). Similar effects were also found for other important interstitials, such as lithium in manganese oxide and carbon in iron. This work establishes a new and improved physical picture for the general diffusion of interstitial atoms in solids and provides new perspectives for us to understand and design diffusion-related material behaviours. For example, the relaxation in metallic glasses 36 was found to correlate with the diffusion of the smallest constituting atoms within some loosely packed/bonded regions. Such diffusion was hypothesized to relate with the vibrations of atomic strings nearby 36 , which is consistent with the current work. It is worth noting that the Brownian motion of nanosized liquid lead (Pb) inclusions within a solid aluminium (Al) was reported 37 . Although the authors found the motion is controlled by the shape of the inclusions and the diffusion of Al along the liquid/solid interface, the agitation force resulting from Al phonons seems to be an important reason for the observed motion.